0.1 Setup

require(devtools)
devtools::install_github("youngser/mbstructure")

0.2 Data Preparation

data(MBconnectome)

0.3 Download data

leftA <- read.table("https://raw.githubusercontent.com/neurodata/graspy/graphmodel/graspy/datasets/drosophila/left_adjacency.csv", header = FALSE, sep = " ")

rightA <- read.table("https://raw.githubusercontent.com/neurodata/graspy/graphmodel/graspy/datasets/drosophila/right_adjacency.csv", header = FALSE, sep = " ")

leftLabels <- read.csv("https://raw.githubusercontent.com/neurodata/graspy/graphmodel/graspy/datasets/drosophila/left_cell_labels.csv", header = FALSE)

rightLabels <- read.csv("https://raw.githubusercontent.com/neurodata/graspy/graphmodel/graspy/datasets/drosophila/right_cell_labels.csv", header = FALSE)

0.4 Embed

X <- as.matrix(rightA)
Xb <- X
Xb[which(Xb > 0)] <- 1

Xase <- irlba(X, nv = 5)

Xl <- Xase$u
Xr <- Xase$v

Xls <- Xase$u / Xase$d
Xrs <- Xase$v / Xase$d

XX <- cbind(Xl, Xr)
XXs <- cbind(Xls, Xrs)

dim(XX)

## [1] 213  10

0.5 Run Urerf

NB: URerF by default normalizes the data before running. This seems to “smooth” the output matrices.

XX.fit  <- Urerf(XX , trees = 500, min.parent = 7, normalizeData = TRUE)
XXs.fit <- Urerf(XXs, trees = 500, min.parent = 7, normalizeData = TRUE)

simMa <- XX.fit$sim
dsimMa <- 1 - simMa

simMb <- XXs.fit$sim
dsimMb <- 1 - simMb

### This is the one we use below
simM <- simMa
dsimM <- dsimMa

simMs <- simMb
dsimMs <- dsimMb

0.6 Plot Affinity Matrices

plot(raster(simM), col = viridis(255), main = "Similarity")

plot(raster(dsimM), col = viridis(255), main = "Dissimilarity")

write.table(file = "similarityMatrix.csv", simM, row.names = FALSE, col.names = FALSE, sep = ",")
write.table(file = "dissimilarityMatrix.csv", 1 - simM, row.names = FALSE, col.names = FALSE, sep = ",")

write.table(file = "similarityMatrixScaledSingularVal.csv", simMs, row.names = FALSE, col.names = FALSE, sep = ",")
write.table(file = "dissimilarityMatrixScaledSingularVal.csv", 1 - simMs, row.names = FALSE, col.names = FALSE, sep = ",")

0.7 Embed the similarity

Y <- irlba(simM, nv = 4)
Ys <- irlba(simMs, nv = 4)

pairs(Y$u, col = scales::alpha(viridis(4)[as.numeric(rightLabels[[1]])], 0.55), pch = 16, cex = 0.75)

pairs(Ys$u, col = scales::alpha(viridis(4)[as.numeric(rightLabels[[1]])], 0.55), pch = 16, cex = 0.75, main = "Scaled by singular values")

1 Mclust

m.fit <-  Mclust(Y$u)
ms.fit <-  Mclust(Ys$u)

plot(m.fit, what = c("BIC"))

plot(ms.fit, what = c("BIC")); title( "Scaled by singular values")

plot(m.fit, what = c("classification"))

plot(ms.fit, what = c("classification")); title( "Scaled by singular values")

mclust::adjustedRandIndex(rightLabels[[1]], m.fit$classification)

## [1] 0.4405913

mclust::adjustedRandIndex(rightLabels[[1]], ms.fit$classification)

## [1] 0.4152494

mi <- 1e4
ARI <- vector("numeric", mi)

for(i in 1:mi){
  set.seed(i)
  ARI[i] <-  
    mclust::adjustedRandIndex(rightLabels[[1]], 
      sample(m.fit$classification, replace= FALSE))
}

plot(histogram(ARI, bin = "Scott", main = "Permutation test of null-distribution"))

2 Hclust

Ydist <- as.dist(dsimM)
h.fit <- hclust(Ydist, method = "average")

plot(h.fit)

plot(as.dendrogram(h.fit))

hclust.labels4 <- cutree(h.fit, k = 4)
hclust.labels5 <- cutree(h.fit, k = 5)
hclust.labels6 <- cutree(h.fit, k = 6)

mclust::adjustedRandIndex(rightLabels[[1]], hclust.labels4)

## [1] 0.6212183

mclust::adjustedRandIndex(rightLabels[[1]], hclust.labels5)

## [1] 0.6425807

mclust::adjustedRandIndex(rightLabels[[1]], hclust.labels6)

## [1] 0.502164

mi <- 1e4
ARIh <- vector("numeric", mi)

for(i in 1:mi){
  set.seed(i)
  ARIh[i] <-  
    mclust::adjustedRandIndex(rightLabels[[1]], 
      sample(hclust.labels5, replace= FALSE))
}

plot(histogram(ARIh, bin = "Scott", main = "Permutation test of null-distribution"))

dend5 <- as.dendrogram(h.fit)

dend5 %>% set("labels", hclust.labels5) %>% 
  color_branches(col = viridis(5)[hclust.labels5]) %>% 
  sort %>%
  plot()

heatmap(as.matrix(Ydist), Rowv = dend5, labCol = hclust.labels5, labRow = hclust.labels5, col = viridis(255))